module SnowpackCompactionAR24Mod

!!! Snowpack compaction process (Anderson 1976) updated by Abolafia-Rosenzweig et al. 2024 (AR24)
!!! Update snow depth via compaction due to destructive metamorphism, overburden, & melt

  use Machine
  use NoahmpVarType
  use ConstantDefineMod

  implicit none

contains

  subroutine SnowpackCompactionAR24(noahmp)

! ------------------------ Code history -----------------------------------
! Original Noah-MP subroutine: COMPACT
! Original code: Guo-Yue Niu and Noah-MP team (Niu et al. 2011)
! Refactered code: C. He, P. Valayamkunnath, & refactor team (He et al. 2023)
! New snow compaction scheme: R. Abolafia-Rosenzweig (Abolafia-Rosenzweig et al. 2024)
! -------------------------------------------------------------------------

    implicit none

    type(noahmp_type), intent(inout) :: noahmp

! local variable
    integer                          :: LoopInd                ! snow layer loop index
    real(kind=kind_noahmp)           :: SnowBurden             ! pressure of overlying snow [kg/m2]
    real(kind=kind_noahmp)           :: SnowCompactAgeExpFac   ! EXPF=exp(-c4*(273.15-TemperatureSoilSnow))
    real(kind=kind_noahmp)           :: TempDiff               ! ConstFreezePoint - TemperatureSoilSnow[K]
    real(kind=kind_noahmp)           :: SnowVoid               ! void (1 - SnowIce - SnowLiqWater)
    real(kind=kind_noahmp)           :: SnowWatTotTmp          ! water mass (ice + liquid) [kg/m2]
    real(kind=kind_noahmp)           :: SnowIceDens            ! partial density of ice [kg/m3]
    real(kind=kind_noahmp)           :: SnowCompactBurdenFac   ! snow overburden compaction parameter [m3/kg]

! --------------------------------------------------------------------
    associate(                                                                        &
              MainTimeStep            => noahmp%config%domain%MainTimeStep           ,& ! in,    noahmp main time step [s]
              TemperatureAirRefHeight => noahmp%forcing%TemperatureAirRefHeight      ,& ! in,    air temperature [K] at reference height
              PressureAirRefHeight    => noahmp%forcing%PressureAirRefHeight         ,& ! in,    air pressure [Pa] at reference height
              TemperatureSoilSnow     => noahmp%energy%state%TemperatureSoilSnow     ,& ! in,    snow and soil layer temperature [K]
              SnowIce                 => noahmp%water%state%SnowIce                  ,& ! in,    snow layer ice [mm]
              SnowLiqWater            => noahmp%water%state%SnowLiqWater             ,& ! in,    snow layer liquid water [mm]
              IndexPhaseChange        => noahmp%water%state%IndexPhaseChange         ,& ! in,    phase change index [0-none;1-melt;2-refreeze]
              SnowIceFracPrev         => noahmp%water%state%SnowIceFracPrev          ,& ! in,    ice fraction in snow layers at previous timestep
              SnowCompactAgingFac1    => noahmp%water%param%SnowCompactAgingFac1     ,& ! in,    snow desctructive metamorphism compaction factor1 [1/s]
              SnowCompactAgingFac2    => noahmp%water%param%SnowCompactAgingFac2     ,& ! in,    snow desctructive metamorphism compaction factor2 [1/k]
              SnowCompactAgingFac3    => noahmp%water%param%SnowCompactAgingFac3     ,& ! in,    snow desctructive metamorphism compaction factor3 
              SnowCompactm            => noahmp%water%param%SnowCompactmAR24         ,& ! in,    snow compaction m parameter for linear sfc temp fitting AR24
              SnowCompactb            => noahmp%water%param%SnowCompactbAR24         ,& ! in,    snow compaction b parameter for linear sfc temp fitting AR24
              SnowCompactPSFC1        => noahmp%water%param%SnowCompactP1AR24        ,& ! in,    lower constraint for SnowCompactBurdenFac for high pressure bin AR24
              SnowCompactPSFC2        => noahmp%water%param%SnowCompactP2AR24        ,& ! in,    lower constraint for SnowCompactBurdenFac for mid pressure bin AR24
              SnowCompactPSFC3        => noahmp%water%param%SnowCompactP3AR24        ,& ! in,    lower constraint for SnowCompactBurdenFac for low pressure bin AR24
              BurdenFacUpper          => noahmp%water%param%BurdenFacUpAR24          ,& ! in,    upper constraint on SnowCompactBurdenFac AR24
              SnowCompactAgingMax     => noahmp%water%param%SnowCompactAgingMax      ,& ! in,    maximum destructive metamorphism compaction [kg/m3]
              SnowViscosityCoeff      => noahmp%water%param%SnowViscosityCoeff       ,& ! in,    snow viscosity coeff [kg s/m2],Anderson1979:0.52e6~1.38e6
              NumSnowLayerNeg         => noahmp%config%domain%NumSnowLayerNeg        ,& ! inout, actual number of snow layers (negative)
              ThicknessSnowSoilLayer  => noahmp%config%domain%ThicknessSnowSoilLayer ,& ! inout, thickness of snow/soil layers [m]
              CompactionSnowAging     => noahmp%water%flux%CompactionSnowAging       ,& ! out,   rate of compaction due to destructive metamorphism [1/s]
              CompactionSnowBurden    => noahmp%water%flux%CompactionSnowBurden      ,& ! out,   rate of compaction of snowpack due to overburden [1/s]
              CompactionSnowMelt      => noahmp%water%flux%CompactionSnowMelt        ,& ! out,   rate of compaction of snowpack due to melt [1/s]
              CompactionSnowTot       => noahmp%water%flux%CompactionSnowTot         ,& ! out,   change in fractional-thickness due to compaction [1/s]
              SnowIceFrac             => noahmp%water%state%SnowIceFrac               & ! out,   fraction of ice in snow layers at current time step
             )
! ----------------------------------------------------------------------

    ! initialization for out-only variables
    CompactionSnowAging(:)  = 0.0
    CompactionSnowBurden(:) = 0.0
    CompactionSnowMelt(:)   = 0.0
    CompactionSnowTot(:)    = 0.0
    SnowIceFrac(:)          = 0.0

    ! start snow compaction
    SnowBurden = 0.0
    ! SnowCompactBurdenFac updated from Abolafia-Rosenzweig et al., 2024
    SnowCompactBurdenFac = SnowCompactm  * TemperatureAirRefHeight + SnowCompactb
    ! pressure-based lower constraints:
    if ( PressureAirRefHeight >= 85000.0 ) then ! high pressure bin
        SnowCompactBurdenFac = max(SnowCompactBurdenFac, SnowCompactPSFC1)
    endif
    if ( (PressureAirRefHeight >= 80000.0) .and. (PressureAirRefHeight < 85000.0) ) then ! mid-pressure bin
        SnowCompactBurdenFac = max(SnowCompactBurdenFac, SnowCompactPSFC2)
    endif
    if ( PressureAirRefHeight < 80000.0 ) then ! low pressure bin
        SnowCompactBurdenFac = max(SnowCompactBurdenFac, SnowCompactPSFC3)
    endif
    ! upper constraint on SnowCompactBurdenFac
    SnowCompactBurdenFac = min(SnowCompactBurdenFac, BurdenFacUpper)

    do LoopInd = NumSnowLayerNeg+1, 0

       SnowWatTotTmp        = SnowIce(LoopInd) + SnowLiqWater(LoopInd)
       SnowIceFrac(LoopInd) = SnowIce(LoopInd) / SnowWatTotTmp
       SnowVoid             = 1.0 - (SnowIce(LoopInd)/ConstDensityIce + SnowLiqWater(LoopInd)/ConstDensityWater) / &
                                    ThicknessSnowSoilLayer(LoopInd)

       ! Allow compaction only for non-saturated node and higher ice lens node.
       if ( (SnowVoid > 0.001) .and. (SnowIce(LoopInd) > 0.1) ) then
          SnowIceDens = SnowIce(LoopInd) / ThicknessSnowSoilLayer(LoopInd)
          TempDiff    = max(0.0, ConstFreezePoint-TemperatureSoilSnow(LoopInd))

          ! Settling/compaction as a result of destructive metamorphism
          SnowCompactAgeExpFac         = exp(-SnowCompactAgingFac2 * TempDiff)
          CompactionSnowAging(LoopInd) = -SnowCompactAgingFac1 * SnowCompactAgeExpFac
          if ( SnowIceDens > SnowCompactAgingMax ) &
             CompactionSnowAging(LoopInd) = CompactionSnowAging(LoopInd) * exp(-46.0e-3*(SnowIceDens-SnowCompactAgingMax))
          if ( SnowLiqWater(LoopInd) > (0.01*ThicknessSnowSoilLayer(LoopInd)) ) &
             CompactionSnowAging(LoopInd) = CompactionSnowAging(LoopInd) * SnowCompactAgingFac3                ! Liquid water term

          ! Compaction due to overburden
          CompactionSnowBurden(LoopInd) = -(SnowBurden + 0.5*SnowWatTotTmp) * &
                                    exp(-0.08*TempDiff-SnowCompactBurdenFac*SnowIceDens) / SnowViscosityCoeff  ! 0.5*SnowWatTotTmp -> self-burden

          ! Compaction occurring during melt
          if ( IndexPhaseChange(LoopInd) == 1 ) then
             CompactionSnowMelt(LoopInd) = max(0.0, (SnowIceFracPrev(LoopInd)-SnowIceFrac(LoopInd)) / &
                                                    max(1.0e-6, SnowIceFracPrev(LoopInd)))
             CompactionSnowMelt(LoopInd) = -CompactionSnowMelt(LoopInd) / MainTimeStep   ! sometimes too large
          else
             CompactionSnowMelt(LoopInd) = 0.0
          endif

          ! Time rate of fractional change in snow thickness (units of s-1)
          CompactionSnowTot(LoopInd) = (CompactionSnowAging(LoopInd) + CompactionSnowBurden(LoopInd) + &
                                        CompactionSnowMelt(LoopInd) ) * MainTimeStep
          CompactionSnowTot(LoopInd) = max(-0.5, CompactionSnowTot(LoopInd))

          ! The change in DZ due to compaction
          ThicknessSnowSoilLayer(LoopInd) = ThicknessSnowSoilLayer(LoopInd) * (1.0 + CompactionSnowTot(LoopInd))
          ThicknessSnowSoilLayer(LoopInd) = max(ThicknessSnowSoilLayer(LoopInd), &
                                            SnowIce(LoopInd)/ConstDensityIce + SnowLiqWater(LoopInd)/ConstDensityWater)

          ! Constrain snow density to a reasonable range (50~500 kg/m3)
          ThicknessSnowSoilLayer(LoopInd) = min( max( ThicknessSnowSoilLayer(LoopInd),&
                                                     (SnowIce(LoopInd)+SnowLiqWater(LoopInd))/500.0 ), &
                                                (SnowIce(LoopInd)+SnowLiqWater(LoopInd))/50.0 )
       endif

       ! Pressure of overlying snow
       SnowBurden = SnowBurden + SnowWatTotTmp

    enddo

    end associate

  end subroutine SnowpackCompactionAR24

end module SnowpackCompactionAR24Mod
